Synthetic data assimilation experiments known as perfect model assimilation are a convenient way to “know” the truth of what is being assimilated. This can help evaluate the data assimilation system and the limits of skill improvements through DA.
Load the rwrfhydro package.
library("rwrfhydro")
## To check rwrfhydro updates run: CheckForUpdates()
# set your own path
rlFile <- '~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03//Route_Link.nc'
# currently there is no "gages" variable
ncdump(rlFile)
## File: ~/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03//Route_Link.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (1):
## linkDim = 403 ;
## variables (24):
## int link(linkDim) ;
## link:long_name = "Link ID" ;
## int from(linkDim) ;
## from:long_name = "From Link ID" ;
## from:coordinates = "time lat lon" ;
## int to(linkDim) ;
## to:long_name = "To Link ID" ;
## to:coordinates = "time lat lon" ;
## float lon(linkDim) ;
## lon:long_name = "longitude of the start node" ;
## lon:units = "degrees_east" ;
## lon:standard_name = "longitude" ;
## float lat(linkDim) ;
## lat:long_name = "latitude of the start node" ;
## lat:units = "degrees_north" ;
## lat:standard_name = "latitude" ;
## float alt(linkDim) ;
## alt:long_name = "Elevation in meters at start node" ;
## alt:standard_name = "height" ;
## alt:units = "m" ;
## alt:positive = "up" ;
## alt:axis = "Z" ;
## int type(linkDim) ;
## type:long_name = "Link type" ;
## type:coordinates = "time lat lon" ;
## int order(linkDim) ;
## order:long_name = "Stream order (Strahler)" ;
## order:coordinates = "time lat lon" ;
## float Qi(linkDim) ;
## Qi:long_name = "Initial flow in link (CMS)" ;
## Qi:coordinates = "time lat lon" ;
## float MusK(linkDim) ;
## MusK:long_name = "Muskingum routing time (s)" ;
## MusK:coordinates = "time lat lon" ;
## float MusX(linkDim) ;
## MusX:long_name = "Muskingum weighting coefficient" ;
## MusX:coordinates = "time lat lon" ;
## float Length(linkDim) ;
## Length:long_name = "Stream length (m)" ;
## Length:coordinates = "time lat lon" ;
## float n(linkDim) ;
## n:long_name = "Manning's roughness" ;
## n:coordinates = "time lat lon" ;
## float So(linkDim) ;
## So:long_name = "Slope (%; drop/length)" ;
## So:coordinates = "time lat lon" ;
## float ChSlp(linkDim) ;
## ChSlp:long_name = "Channel side slope (%; drop/length)" ;
## ChSlp:coordinates = "time lat lon" ;
## float BtmWdth(linkDim) ;
## BtmWdth:long_name = "Bottom width of channel" ;
## BtmWdth:coordinates = "time lat lon" ;
## float LkHZArea(linkDim) ;
## LkHZArea:long_name = "Maximum lake area (sq. m)" ;
## LkHZArea:coordinates = "time lat lon" ;
## float LkMxH(linkDim) ;
## LkMxH:long_name = "Maximum lake elevation (m)" ;
## LkMxH:coordinates = "time lat lon" ;
## float WeirC(linkDim) ;
## WeirC:long_name = "Weir coefficient" ;
## WeirC:coordinates = "time lat lon" ;
## float WeirL(linkDim) ;
## WeirL:long_name = "Weir length (m)" ;
## WeirL:coordinates = "time lat lon" ;
## float OrificeC(linkDim) ;
## OrificeC:long_name = "Orifice coefficient" ;
## OrificeC:coordinates = "time lat lon" ;
## float OrificeA(linkDim) ;
## OrificeA:long_name = "Orifice croess-sectional area (sq. m)" ;
## OrificeA:coordinates = "time lat lon" ;
## float OrificeE(linkDim) ;
## OrificeE:long_name = "Orifice elevation (m)" ;
## OrificeE:coordinates = "time lat lon" ;
## float time() ;
## time:standard_name = "time" ;
## time:long_name = "time of measurement" ;
## time:units = "days since 2000-01-01 00:00:00" ;
##
## // global attributes (2):
## :featureType = "point"
## :history = "Created Thu Sep 03 09:30:05 2015"
# This function returns a function which plots the data
PlotRouteLink <- VisualizeRouteLink(rlFile)
PlotRouteLink() # basic plot
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.940325,-105.430153&zoom=8&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
args(PlotRouteLink) # what are the options? unfortunately I havent really documented these.
## function (location = c(lon = mean(range(rl$lon)), lat = mean(range(rl$lat))),
## indices = FALSE, comIds = FALSE, zoom = if (length(gageZoom)) 13 else 8,
## source = "google", maptype = "terrain", padPlot = 0.01, linkColor = "blue",
## textColor = "darkred", gageColor = "cyan", gageShape = 9,
## gageZoom = NULL, plotPath = NULL, plotType = "pdf", doPlot = TRUE)
## NULL
PlotRouteLink(zoom=10, comIds=TRUE) # links on the edge of the domain arent fully plotted... because there's no "to" lat/lon.
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=39.940325,-105.430153&zoom=10&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
Define the synthetic gages to be used for the domain. In this example put gages on ALL the links.
rlLink <- ncdump(rlFile, 'link', quiet=TRUE)
gageIds <- paste0('g',formatC(rlLink,width=3,flag='0'))
newCopyId <- 'gagesAllLinks'
newFile <- AddRouteLinkGage(rlFile, gageIds, rlLink, new=newCopyId, overwrite=TRUE)
print(newFile)
## [1] "/Users/jamesmcc/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc"
ncdump(newFile)
## File: /Users/jamesmcc/WRF_Hydro/DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (2):
## linkDim = 403 ;
## IDLength = 15 ;
## variables (26):
## integer IDLength(IDLength) ;
## IDLength:long_name = "IDLength"
## int link(linkDim) ;
## link:long_name = "Link ID" ;
## int from(linkDim) ;
## from:long_name = "From Link ID" ;
## from:coordinates = "time lat lon" ;
## int to(linkDim) ;
## to:long_name = "To Link ID" ;
## to:coordinates = "time lat lon" ;
## float lon(linkDim) ;
## lon:long_name = "longitude of the start node" ;
## lon:units = "degrees_east" ;
## lon:standard_name = "longitude" ;
## float lat(linkDim) ;
## lat:long_name = "latitude of the start node" ;
## lat:units = "degrees_north" ;
## lat:standard_name = "latitude" ;
## float alt(linkDim) ;
## alt:long_name = "Elevation in meters at start node" ;
## alt:standard_name = "height" ;
## alt:units = "m" ;
## alt:positive = "up" ;
## alt:axis = "Z" ;
## int type(linkDim) ;
## type:long_name = "Link type" ;
## type:coordinates = "time lat lon" ;
## int order(linkDim) ;
## order:long_name = "Stream order (Strahler)" ;
## order:coordinates = "time lat lon" ;
## float Qi(linkDim) ;
## Qi:long_name = "Initial flow in link (CMS)" ;
## Qi:coordinates = "time lat lon" ;
## float MusK(linkDim) ;
## MusK:long_name = "Muskingum routing time (s)" ;
## MusK:coordinates = "time lat lon" ;
## float MusX(linkDim) ;
## MusX:long_name = "Muskingum weighting coefficient" ;
## MusX:coordinates = "time lat lon" ;
## float Length(linkDim) ;
## Length:long_name = "Stream length (m)" ;
## Length:coordinates = "time lat lon" ;
## float n(linkDim) ;
## n:long_name = "Manning's roughness" ;
## n:coordinates = "time lat lon" ;
## float So(linkDim) ;
## So:long_name = "Slope (%; drop/length)" ;
## So:coordinates = "time lat lon" ;
## float ChSlp(linkDim) ;
## ChSlp:long_name = "Channel side slope (%; drop/length)" ;
## ChSlp:coordinates = "time lat lon" ;
## float BtmWdth(linkDim) ;
## BtmWdth:long_name = "Bottom width of channel" ;
## BtmWdth:coordinates = "time lat lon" ;
## float LkHZArea(linkDim) ;
## LkHZArea:long_name = "Maximum lake area (sq. m)" ;
## LkHZArea:coordinates = "time lat lon" ;
## float LkMxH(linkDim) ;
## LkMxH:long_name = "Maximum lake elevation (m)" ;
## LkMxH:coordinates = "time lat lon" ;
## float WeirC(linkDim) ;
## WeirC:long_name = "Weir coefficient" ;
## WeirC:coordinates = "time lat lon" ;
## float WeirL(linkDim) ;
## WeirL:long_name = "Weir length (m)" ;
## WeirL:coordinates = "time lat lon" ;
## float OrificeC(linkDim) ;
## OrificeC:long_name = "Orifice coefficient" ;
## OrificeC:coordinates = "time lat lon" ;
## float OrificeA(linkDim) ;
## OrificeA:long_name = "Orifice croess-sectional area (sq. m)" ;
## OrificeA:coordinates = "time lat lon" ;
## float OrificeE(linkDim) ;
## OrificeE:long_name = "Orifice elevation (m)" ;
## OrificeE:coordinates = "time lat lon" ;
## float time() ;
## time:standard_name = "time" ;
## time:long_name = "time of measurement" ;
## time:units = "days since 2000-01-01 00:00:00" ;
## char gages(linkDim,IDLength) ;
## gages:units = "usgsId" ;
## gages:_FillValue = " " ;
##
## // global attributes (2):
## :featureType = "point"
## :history = "Created Thu Sep 03 09:30:05 2015"
print(ncdump(newFile, 'gages', quiet=TRUE))
## [1] " g255" " g266" " g213"
## [4] " g253" " g277" " g208"
## [7] " g259" " g240" " g241"
## [10] " g257" " g217" " g216"
## [13] " g238" " g228" " g247"
## [16] " g229" " g231" " g258"
## [19] " g069" " g230" " g061"
## [22] " g244" " g063" " g222"
## [25] " g263" " g262" " g066"
## [28] " g256" " g252" " g115"
## [31] " g141" " g140" " g225"
## [34] " g103" " g067" " g065"
## [37] " g091" " g090" " g260"
## [40] " g187" " g227" " g146"
## [43] " g112" " g210" " g264"
## [46] " g133" " g132" " g163"
## [49] " g226" " g128" " g221"
## [52] " g092" " g215" " g110"
## [55] " g177" " g142" " g139"
## [58] " g234" " g199" " g164"
## [61] " g106" " g123" " g169"
## [64] " g118" " g209" " g170"
## [67] " g105" " g200" " g197"
## [70] " g126" " g223" " g186"
## [73] " g153" " g120" " g147"
## [76] " g175" " g107" " g179"
## [79] " g224" " g096" " g124"
## [82] " g154" " g089" " g151"
## [85] " g211" " g176" " g203"
## [88] " g138" " g167" " g099"
## [91] " g125" " g155" " g152"
## [94] " g150" " g181" " g207"
## [97] " g161" " g097" " g190"
## [100] " g058" " g280" " g278"
## [103] " g050" " g144" " g174"
## [106] " g098" " g188" " g059"
## [109] " g282" " g281" " g055"
## [112] " g113" " g204" " g104"
## [115] " g166" " g165" " g189"
## [118] " g121" " g056" " g052"
## [121] " g178" " g172" " g171"
## [124] " g300" " g295" " g101"
## [127] " g191" " g094" " g093"
## [130] " g316" " g122" " g057"
## [133] " g180" " g049" " g303"
## [136] " g302" " g325" " g095"
## [139] " g317" " g314" " g343"
## [142] " g053" " g212" " g040"
## [145] " g341" " g100" " g194"
## [148] " g087" " g082" " g332"
## [151] " g330" " g198" " g195"
## [154] " g319" " g086" " g083"
## [157] " g235" " g201" " g326"
## [160] " g318" " g085" " g205"
## [163] " g108" " g196" " g320"
## [166] " g088" " g310" " g116"
## [169] " g145" " g304" " g109"
## [172] " g074" " g232" " g068"
## [175] " g192" " g184" " g183"
## [178] " g084" " g305" " g296"
## [181] " g071" " g102" " g070"
## [184] " g292" " g288" " g315"
## [187] " g313" " g290" " g182"
## [190] " g080" " g078" " g173"
## [193] " g299" " g297" " g168"
## [196] " g162" " g289" " g157"
## [199] " g156" " g079" " g309"
## [202] " g308" " g143" " g075"
## [205] " g298" " g136" " g294"
## [208] " g285" " g248" " g119"
## [211] " g117" " g076" " g137"
## [214] " g261" " g291" " g322"
## [217] " g321" " g287" " g283"
## [220] " g249" " g329" " g242"
## [223] " g336" " g111" " g073"
## [226] " g038" " g324" " g323"
## [229] " g032" " g286" " g251"
## [232] " g345" " g254" " g340"
## [235] " g051" " g338" " g337"
## [238] " g062" " g339" " g328"
## [241] " g327" " g293" " g352"
## [244] " g348" " g347" " g346"
## [247] " g361" " g054" " g021"
## [250] " g243" " g018" " g047"
## [253] " g334" " g333" " g265"
## [256] " g335" " g359" " g006"
## [259] " g356" " g355" " g033"
## [262] " g127" " g378" " g376"
## [265] " g369" " g239" " g366"
## [268] " g365" " g362" " g360"
## [271] " g358" " g357" " g349"
## [274] " g214" " g077" " g012"
## [277] " g311" " g048" " g301"
## [280] " g039" " g034" " g030"
## [283] " g025" " g276" " g019"
## [286] " g273" " g268" " g267"
## [289] " g131" " g379" " g377"
## [292] " g375" " g374" " g245"
## [295] " g371" " g114" " g368"
## [298] " g367" " g237" " g364"
## [301] " g363" " g354" " g353"
## [304] " g350" " g284" " g014"
## [307] " g081" " g193" " g031"
## [310] " g312" " g046" " g041"
## [313] " g036" " g035" " g159"
## [316] " g028" " g027" " g023"
## [319] " g022" " g020" " g275"
## [322] " g015" " g270" " g269"
## [325] " g388" " g382" " g381"
## [328] " g380" " g373" " g372"
## [331] " g236" " g233" " g016"
## [334] " g351" " g220" " g218"
## [337] " g342" " g386" " g064"
## [340] " g060" " g307" " g306"
## [343] " g045" " g042" " g037"
## [346] " g160" " g389" " g029"
## [349] " g026" " g149" " g274"
## [352] " g017" " g272" " g013"
## [355] " g008" " g135" " g005"
## [358] " g387" " g130" " g385"
## [361] " g384" " g383" " g250"
## [364] " g246" " g370" " g391"
## [367] " g219" " g344" " g399"
## [370] " g206" " g331" " g202"
## [373] " g072" " g395" " g394"
## [376] " g185" " g393" " g044"
## [379] " g043" " g129" " g390"
## [382] " g158" " g024" " g279"
## [385] " g148" " g403" " g402"
## [388] " g401" " g400" " g271"
## [391] " g398" " g397" " g396"
## [394] " g011" " g010" " g009"
## [397] " g392" " g007" " g134"
## [400] " g004" " g003" " g002"
## [403] " g001"
Set the symlink to the new Route_Link.nc file.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek/DOMAIN> l
lrwxrwxrwx 2 jamesmcc rap 93 Sep 9 09:42 Route_Link.nc ->
../../DOMAIN_library/Boulder_Creek_100m_1km_2sqkm_full_2015_09_03/Route_Link.gagesAllLinks.nc
Set up several model parameters in the namelists.
namelist.hrldas:
START_YEAR = 2012
START_MONTH = 08
START_DAY = 01
START_HOUR = 00
START_MIN = 00
KHOUR = 48
FORCING_TIMESTEP = 3600
NOAH_TIMESTEP = 900
OUTPUT_TIMESTEP = 7200
FORC_TYP = 4
hydro.namelist:
out_dt = 15 ! minutes
DTRT_CH = 60
DTRT_TER = 6
SUBRTSWCRT = 1
OVRTSWCRT = 1
CHANRTSWCRT = 1
channel_option = 2
route_link_f = "DOMAIN/Route_Link.nc"
Run without nudging and using precip doubling. The wrf_hydro.noNudging_doublePrecip.exe binary was built with the following environment variables: PRECIP_DOUBLE=1 WRF_HYDRO_NUDGING=0 using the daBranch of wrf_hydro_model.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.noNudging_doublePrecip.exe
. . .
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmo
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmo
Check chanobs and frxst pts dimensions
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek/run.pmo> wc -l frxst_pts_out.txt
77779 frxst_pts_out.txt ## = 403stns*(48hours*4output/hr+1output@0)
frxst <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmo/frxst_pts_out.txt", stId='character')
ggplot2::ggplot(subset(frxst, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')), ## all the action is early on
ggplot2::aes(x=POSIXct, y=q_cms, color=st_id)) +
ggplot2::geom_line() +
ggplot2::scale_color_discrete(guide='none') +
ggplot2::theme_bw()
Plot it by order
rlOrder <- ncdump(rlFile, 'order', quiet=TRUE)
names(rlOrder) <- paste0('g',formatC(rlLink,width=3,flag='0'))
frxst$order <- rlOrder[trimws(frxst$st_id)]
ggplot2::ggplot(subset(frxst, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')),
ggplot2::aes(x=POSIXct, y=q_cms, color=st_id)) +
ggplot2::geom_line() +
ggplot2::scale_color_discrete(guide='none') +
ggplot2::facet_wrap(~order, ncol=1, scales='free_y') +
ggplot2::theme_bw()
chanObsFiles <- list.files('~/WRF_Hydro/Col_Bldr_Creek/run.pmo/', pattern = 'CHANOBS_DOMAIN', full.names=TRUE)
sliceResolutionMin <- 15
outputDir <- '~/WRF_Hydro/Col_Bldr_Creek/run.pmo/timeSlicePMO/'
dir.create(outputDir)
sliceFiles <- ChanObsToTimeSlice(chanObsFiles, sliceResolutionMin, outputDir)
Put these in place
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> rm nudgingTimeSliceObs
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ln -sf run.pmo/timeSlicePMO/ nudgingTimeSliceObs
min(ncdump(rlFile, 'Length', quiet=TRUE)) ## use R less than this /2.
## [1] 50
gageParams <- data.frame( gageId = paste0('g',formatC(1:403,width=3,flag='0')), stringsAsFactors=FALSE)
gageParams$R=24
gageParams$G=1
gageParams$tau=20
MkNudgingParams(gageId=gageParams$gageId, R=gageParams$R,
G=gageParams$G, tau=gageParams$tau,
outFile='~/WRF_Hydro/Col_Bldr_Creek/run.pmo/nudgingParams.PMOallGages.nc',
overwrite=TRUE)
## File: ~/WRF_Hydro/Col_Bldr_Creek/run.pmo/nudgingParams.PMOallGages.nc
## ( NC_FORMAT_NETCDF4 ):
## dimensions (2):
## stationIdStrLen = 15 ;
## stationIdInd = UNLIMITED ; // (403 currently)
## variables (4):
## char stationId(stationIdInd,stationIdStrLen) (Chunking: [15,1]) ;
## stationId:units = "-" ;
## stationId:long_name = "USGS station identifer" ;
## float R(stationIdInd) (Chunking: [1]) ;
## R:units = "meters" ;
## R:long_name = "Radius of influence in meters" ;
## float G(stationIdInd) (Chunking: [1]) ;
## G:units = "-" ;
## G:long_name = "Amplitude of nudging" ;
## float tau(stationIdInd) (Chunking: [1]) ;
## tau:units = "minutes" ;
## tau:long_name = "Time tapering parameter half window size in minutes" ;
Then put the param file in place.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ln -s run.pmo/nudgingParams.PMOallGages.nc nudgingParams.nc
The wrf_hydro.nudging.exe binary was built with the following environment variables: PRECIP_DOUBLE=0 and WRF_HYDRO_NUDGING=1 using the daBranch of wrf_hydro_model.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.nudging.exe
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmoNudge
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmoNudge
The wrf_hydro.noNudging.exe binary was built with the following environment variables: PRECIP_DOUBLE=0 and WRF_HYDRO_NUDGING=0 using the daBranch of wrf_hydro_model.
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mpirun -np 1 ./wrf_hydro.noNudging.exe
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> mkdir run.pmoNoNudge
jamesmcc@hydro-c1:~/WRF_Hydro/Col_Bldr_Creek> ./cleanup.sh run.pmoNoNudge
77779 frxst_pts_out.txt ## = 403stns*(48hours*4output/hr+1output@0)
frxstNoNudge <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmoNoNudge/frxst_pts_out.txt", stId='character')
frxstNudge <- ReadFrxstPts("~/WRF_Hydro/Col_Bldr_Creek/run.pmoNudge/frxst_pts_out.txt", stId='character')
## make sure the abscissae columns are equal before taking the differences
if(all(plyr::laply(NamedList(names(frxst)), function(colName) all(frxst[[colName]]==frxstNoNudge[[colName]]))[-6:-7]))
{
forcErr <- frxst
forcErr$`q err obs-model (cms)` <- frxst$q_cms - frxstNoNudge$q_cms
forcErr$`q err obs-model (cfs)` <- frxst$q_cfs - frxstNoNudge$q_cfs
forcErr$q_cms <- forcErr$q_cfs <- NULL
forcErr$error <- 'forcing'
} else warning('Abscissae do not line up with frxstNoNudge, please investigate', immediate.=TRUE)
if(all(plyr::laply(NamedList(names(frxst)), function(colName) all(frxst[[colName]]==frxstNudge[[colName]]))[-6:-7]))
{
nudgeErr <- frxst
nudgeErr$`q err obs-model (cms)` <- frxst$q_cms - frxstNudge$q_cms
nudgeErr$`q err obs-model (cfs)` <- frxst$q_cfs - frxstNudge$q_cfs
nudgeErr$q_cms <- nudgeErr$q_cfs <- NULL
nudgeErr$error <- 'nudging'
} else warning('Abscissae do not line up with frxstNudge, please investigate', immediate.=TRUE)
err <- rbind(forcErr, nudgeErr)
ggplot2::ggplot(subset(err, POSIXct < as.POSIXct('2012-08-01 06:00:00', tz='UTC')), ## all the action is early on
ggplot2::aes(x=POSIXct, y=`q err obs-model (cms)`, color=st_id)) +
ggplot2::geom_line() +
ggplot2::scale_color_discrete(guide='none') +
ggplot2::facet_wrap(~error, ncol=1) +
ggplot2::theme_bw()